##### LAGGED DEPENDENT VARIABLE #### 

library(tidyverse)
library(broom)
library(fixest)
library(ggthemes)

source("code/01-importdata.R")
source("code/02-mainresults.R")


## Re-run models with lagged dependent variable, remove country-specific time trends

ldv_laws <- feols(c(scaled_climate_laws)
                    ~ l(x = scaled_climate_laws, 1) + annual_W_lag1 +
                      sw(annual_average_temp_lag1, annual_average_temp_lag2,
                                                    annual_average_temp_lag3,
                                                    hottest_season_temp_lag1, hottest_season_temp_lag2,
                                                    hottest_season_temp_lag3) |
                      # iso3c[year] + year, # trends
                    iso3c + year, # twoway
                    panel.id = ~iso3c+year,
                    data = analysis_df)

ldv_ecp <- feols(c(scaled_ecp_zeros)
                  ~ l(x = scaled_ecp_zeros, 1) + annual_W_lag1 +
                    sw(annual_average_temp_lag1, annual_average_temp_lag2,
                       annual_average_temp_lag3,
                       hottest_season_temp_lag1, hottest_season_temp_lag2,
                       hottest_season_temp_lag3) |
                   # iso3c[year] + year, # trends
                   iso3c + year, # twoway
                  panel.id = ~iso3c+year,
                  data = analysis_df)

ldv_fit <- feols(c(scaled_fit)
                 ~ l(x = scaled_fit, 1) + annual_W_lag1 +
                   sw(annual_average_temp_lag1, annual_average_temp_lag2,
                      annual_average_temp_lag3,
                      hottest_season_temp_lag1, hottest_season_temp_lag2,
                      hottest_season_temp_lag3) |
                   # iso3c[year] + year, # trends
                   iso3c + year, # twoway
                 panel.id = ~iso3c+year,
                 data = analysis_df)

ldv_quota <- feols(c(scaled_quota)
                 ~ l(x = scaled_quota, 1) + annual_W_lag1 +
                   sw(annual_average_temp_lag1, annual_average_temp_lag2,
                      annual_average_temp_lag3,
                      hottest_season_temp_lag1, hottest_season_temp_lag2,
                      hottest_season_temp_lag3) |
                   # iso3c[year] + year, # trends
                   iso3c + year, # twoway
                 panel.id = ~iso3c+year,
                 data = analysis_df)

ldv_renew <- feols(c(scaled_elec_renew)
                 ~ l(x = scaled_elec_renew, 1) + annual_W_lag1 +
                   sw(annual_average_temp_lag1, annual_average_temp_lag2,
                      annual_average_temp_lag3,
                      hottest_season_temp_lag1, hottest_season_temp_lag2,
                      hottest_season_temp_lag3) |
                   # iso3c[year] + year, # trends
                   iso3c + year, # twoway
                 panel.id = ~iso3c+year,
                 data = analysis_df)

ldv_delegates <- feols(c(scaled_delegates)
                   ~ l(x = scaled_delegates, 1) + annual_W_lag1 +
                     sw(annual_average_temp_lag1, annual_average_temp_lag2,
                        annual_average_temp_lag3,
                        hottest_season_temp_lag1, hottest_season_temp_lag2,
                        hottest_season_temp_lag3) |
                     # iso3c[year] + year, # trends
                     iso3c + year, # twoway
                   panel.id = ~iso3c+year,
                   data = analysis_df)

ldv_gcg <- feols(c(scaled_gcg)
                       ~ l(x = scaled_gcg, 1) + annual_W_lag1 +
                         sw(annual_average_temp_lag1, annual_average_temp_lag2,
                            annual_average_temp_lag3,
                            hottest_season_temp_lag1, hottest_season_temp_lag2,
                            hottest_season_temp_lag3) |
                   # iso3c[year] + year, # trends
                   iso3c + year, # twoway
                       panel.id = ~iso3c+year,
                       data = analysis_df)

ldv_total_provided <- feols(c(scaled_cf_total_provided)
                 ~ l(x = scaled_cf_total_provided, 1) + annual_W_lag1 +
                   sw(annual_average_temp_lag1, annual_average_temp_lag2,
                      annual_average_temp_lag3,
                      hottest_season_temp_lag1, hottest_season_temp_lag2,
                      hottest_season_temp_lag3) |
                   # iso3c[year] + year, # trends
                   iso3c + year, # twoway
                 panel.id = ~iso3c+year,
                 data = analysis_df)


ldv_prin_provided <- feols(c(scaled_cf_principal_provided)
                 ~ l(x = scaled_cf_principal_provided, 1) + annual_W_lag1 +
                   sw(annual_average_temp_lag1, annual_average_temp_lag2,
                      annual_average_temp_lag3,
                      hottest_season_temp_lag1, hottest_season_temp_lag2,
                      hottest_season_temp_lag3) |
                   # iso3c[year] + year, # trends
                   iso3c + year, # twoway
                 panel.id = ~iso3c+year,
                 data = analysis_df)


ldv_total_received <- feols(c(scaled_cf_total_received)
                 ~ l(x = scaled_cf_total_received, 1) + annual_W_lag1 +
                   sw(annual_average_temp_lag1, annual_average_temp_lag2,
                      annual_average_temp_lag3,
                      hottest_season_temp_lag1, hottest_season_temp_lag2,
                      hottest_season_temp_lag3) |
                   # iso3c[year] + year, # trends
                   iso3c + year, # twoway
                 panel.id = ~iso3c+year,
                 data = analysis_df)


ldv_prin_received <- feols(c(scaled_cf_principal_received)
                 ~ l(x = scaled_cf_principal_received, 1) + annual_W_lag1 +
                   sw(annual_average_temp_lag1, annual_average_temp_lag2,
                      annual_average_temp_lag3,
                      hottest_season_temp_lag1, hottest_season_temp_lag2,
                      hottest_season_temp_lag3) |
                   # iso3c[year] + year, # trends
                   iso3c + year, # twoway
                 panel.id = ~iso3c+year,
                 data = analysis_df)

# Combine the results to be in the same format as the main results
ldv_list <- c(as.list(ldv_laws),
              as.list(ldv_ecp),
              as.list(ldv_renew),
              as.list(ldv_fit),
              as.list(ldv_quota),
              as.list(ldv_delegates),
              as.list(ldv_gcg),
              as.list(ldv_total_provided),
              as.list(ldv_prin_provided),
              as.list(ldv_total_received),
              as.list(ldv_prin_received))

ldv_results <-  setNames(as.data.frame(matrix(NA, nrow = 66, ncol = 10)), 
                            c("outcome", "treatment", "beta", "std_error",
                              "t_stat", "p_value", "nobs", "n_countries",
                              "min_year", "max_year"))

for (i in 1:66){
  ldv_results[i, 1] <- ldv_list[[i]]$fml %>% gsub("\\~", "", .) %>% .[2] # outcome
  ldv_results[i, 2] <- ldv_list[[i]]$fml %>% gsub("\\~", "", .) %>% .[3] # treatment
  ldv_results[i, 3] <- ldv_list[[i]] %>% tidy() %>% slice(3) %>% filter(!grepl("_W_",term)) %>% select(2) # beta
  ldv_results[i, 4] <- ldv_list[[i]] %>% tidy() %>% slice(3) %>% filter(!grepl("_W_",term)) %>% select(3) # std error
  ldv_results[i, 5] <- ldv_list[[i]] %>% tidy() %>% slice(3) %>% filter(!grepl("_W_",term)) %>% select(4) # t-stat
  ldv_results[i, 6] <- ldv_list[[i]] %>% tidy() %>% slice(3) %>% filter(!grepl("_W_",term)) %>% select(5) # p-val
  ldv_results[i, 7] <- ldv_list[[i]]$nobs # number of observations
  ldv_results[i, 8] <- length(unique(ldv_list[[i]]$fixef_id$iso3c)) # number of countries
  ldv_results[i, 9] <- min(attr(ldv_list[[i]]$fixef_id$year, "fixef_names")) # start year of analysis
  ldv_results[i, 10] <- max(attr(ldv_list[[i]]$fixef_id$year, "fixef_names")) # end year of analysis
}

ldv_results$treatment <- str_sub(ldv_results$treatment, start = -24L, -1L)

## Combine results from baseline trends and LDV models

ldv_trends_combine <- cbind.data.frame(ldv_results %>% 
            select(outcome, treatment, ldv_t_stat = t_stat),
          results_trends %>% 
            select(baseline_t_stat = t_stat))

ldv_trends_combine %>% 
  ggplot(., aes(baseline_t_stat, ldv_t_stat)) +
  coord_cartesian(xlim = c(-5,5), ylim = c(-5,5)) +
  geom_abline(slope = 1, intercept = 0, linetype = 2, color = "grey50") +
  geom_vline(xintercept = c(-2.84, -1.96, 1.96, 2.84),
             color = c("#79b321", "#456613", "#456613", "#79b321"),
             linetype = 2,
             size = c(1.5, 1, 1, 1.5)) +
  geom_hline(yintercept = c(-2.84, -1.96, 1.96, 2.84),
             color = c("#79b321", "#456613", "#456613", "#79b321"),
             linetype = 2,
             size = c(1.5, 1, 1, 1.5)) +
  geom_point(size = 2) +
  scale_x_continuous(breaks = c(-5, -2.84, -1.96, -1, 0, 1, 1.96, 2.84, 5)) +
  scale_y_continuous(breaks = c(-5, -2.84, -1.96, -1, 0, 1, 1.96, 2.84, 5)) +
  labs(x="T-statistics in baseline model", 
       y="T-statistics with lagged DV") +
  theme_tufte(base_family = "Helvetica") +
  theme(legend.position = "none",
        panel.background = element_rect(fill=NA),
        # panel.grid.major.x = element_line(colour = "grey90"),
        # panel.grid.minor.x = element_line(colour = "grey90"),
        # panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=16))

ggsave("text/robustness_ldv.pdf", units = "in", height = 8, width = 10)

